#include <admodel.h>

  extern "C"  {
    void ad_boundf(int i);
  }
#include <msm12.htp>

model_data::model_data(int argc,char * argv[]) : ad_comm(argc,argv)
{
MIXM2=.9;
niter=4;
  mode.allocate("mode");
  styr.allocate("styr");
  nyrs.allocate("nyrs");
  nspp.allocate("nspp");
  nages.allocate("nages");
  obs_catch.allocate(1,nspp,1,nyrs,1,nages,"obs_catch");
  wt.allocate(1,nspp,1,nspp,1,nages,1,nages,"wt");
  food_in.allocate(1,nspp,1,nyrs,1,nages,"food_in");
  avg_growth.allocate(1,nspp,1,nyrs,1,nages,"avg_growth");
  other_food.allocate(1,2,"other_food");
  suit_main.allocate(1,nspp,1,nspp,1,nages,1,nages,"suit_main");
  suit_other.allocate(1,2,1,nages,"suit_other");
  M1.allocate(1,nspp,1,nages,"M1");
  suit_transp.allocate(1,nspp,1,nspp,1,nages,1,nages,"suit_transp");
  pred_food_ration.allocate(1,nspp,1,nyrs,1,nages);
  suit_wt.allocate(1,nspp,1,nages);
  trend.allocate(1,nyrs);
 trend.fill_seqadd(-1,2/(nyrs-1));
  p_obs.allocate(1,nspp,1,nyrs,1,nages);
  tc_obs.allocate(1,nspp,1,nyrs);
   tau = 50.;
  
   for (k=1;k<=nspp;k++)
   {
     for (i=1;i<=nyrs;i++)
     {
       tc_obs(k,i)= sum(obs_catch(k,i));
       p_obs(k,i) = obs_catch(k,i)/(tc_obs(k,i)+.01);
       
     }
   }
  n_srv_biom.allocate(1,nspp,"n_srv_biom");
  yrs_srv_biom.allocate(1,nspp,1,n_srv_biom,"yrs_srv_biom");
  srv_biom.allocate(1,nspp,1,n_srv_biom,"srv_biom");
  srv_biom_se.allocate(1,nspp,1,n_srv_biom,"srv_biom_se");
  srv_biom_lse.allocate(1,nspp,1,n_srv_biom);
 srv_biom_lse = elem_div(srv_biom_se,srv_biom);
 srv_biom_lse = sqrt(log(square(srv_biom_lse) + 1.));
 cout<<srv_biom_lse<<endl;
  n_srv_age.allocate(1,nspp,"n_srv_age");
  srv_age_type.allocate(1,nspp,"srv_age_type");
  yrs_srv_age.allocate(1,nspp,1,n_srv_age,"yrs_srv_age");
  srv_age_bins.allocate(1,nspp,"srv_age_bins");
  srv_age_n.allocate(1,nspp,1,n_srv_age,"srv_age_n");
  srv_age.allocate(1,nspp,1,n_srv_age,1,srv_age_bins,"srv_age");
for (k=1;k<=nspp;k++) for (i=1;i<=n_srv_age(k);i++) srv_age(k,i) /= sum(srv_age(k,i));
  srv_age_err.allocate(1,nspp,1,nages,1,srv_age_bins,"srv_age_err");
  test_read.allocate("test_read");
 if (test_read != 12345) {cout<<"Read file error"<<endl<<test_read<<endl<<srv_age_err<<endl;exit(1);}
}

void model_parameters::initializationfunction(void)
{
  srv_sel_inf.set_initial_value(4.5);
  srv_sel_slp.set_initial_value(.9);
  mean_rec.set_initial_value(2.);
  log_srv_q.set_initial_value(0.);
  mean_F.set_initial_value(.2);
  sel_coff.set_initial_value(1.);
  alpha.set_initial_value(1);
}

model_parameters::model_parameters(int sz,int argc,char * argv[]) : 
 model_data(argc,argv) , function_minimizer(sz)
{
  initializationfunction();
  mean_rec.allocate(1,nspp,1,"mean_rec");
  mean_F.allocate(1,nspp,.001,2,1,"mean_F");
  rec_dev.allocate(1,nspp,1,nyrs,-10,10,2,"rec_dev");
  init_dev.allocate(1,nspp,2,nages,-10,10,2,"init_dev");
  F_dev.allocate(1,nspp,1,nyrs,-4,4,2,"F_dev");
  alpha.allocate(.49,1.01,-3,"alpha");
 nselages=8;
  sel_coff.allocate(1,nspp,1,nselages,4,"sel_coff");
  sel.allocate(1,nspp,1,nages,"sel");
  #ifndef NO_AD_INITIALIZE
    sel.initialize();
  #endif
  M2_pen.allocate("M2_pen");
  #ifndef NO_AD_INITIALIZE
  M2_pen.initialize();
  #endif
  R.allocate(1,nspp,1,nyrs,"R");
  #ifndef NO_AD_INITIALIZE
    R.initialize();
  #endif
  biomass.allocate(1,nspp,1,nyrs,"biomass");
  #ifndef NO_AD_INITIALIZE
    biomass.initialize();
  #endif
  depletion.allocate(1,nspp,"depletion");
  p_hat.allocate(1,nspp,1,nyrs,1,nages,"p_hat");
  #ifndef NO_AD_INITIALIZE
    p_hat.initialize();
  #endif
  tc_hat.allocate(1,nspp,1,nyrs,"tc_hat");
  #ifndef NO_AD_INITIALIZE
    tc_hat.initialize();
  #endif
  log_srv_q.allocate(1,nspp,3,"log_srv_q");
  srv_sel_inf.allocate(1,nspp,-3,"srv_sel_inf");
  srv_sel_slp.allocate(1,nspp,-5,"srv_sel_slp");
  srv_sel_coff.allocate(1,nspp,1,nselages,3,"srv_sel_coff");
  avg_srv_sel.allocate(1,nspp,"avg_srv_sel");
  #ifndef NO_AD_INITIALIZE
    avg_srv_sel.initialize();
  #endif
  avg_sel.allocate(1,nspp,"avg_sel");
  #ifndef NO_AD_INITIALIZE
    avg_sel.initialize();
  #endif
  srv_q.allocate(1,nspp,"srv_q");
  #ifndef NO_AD_INITIALIZE
    srv_q.initialize();
  #endif
  srv_biom_hat.allocate(1,nspp,1,n_srv_biom,"srv_biom_hat");
  #ifndef NO_AD_INITIALIZE
    srv_biom_hat.initialize();
  #endif
  srv_sel.allocate(1,nspp,1,nages,"srv_sel");
  #ifndef NO_AD_INITIALIZE
    srv_sel.initialize();
  #endif
  srv_age_hat.allocate(1,nspp,1,n_srv_age,1,srv_age_bins,"srv_age_hat");
  #ifndef NO_AD_INITIALIZE
    srv_age_hat.initialize();
  #endif
  catch_hat.allocate(1,nspp,1,nyrs,1,nages,"catch_hat");
  #ifndef NO_AD_INITIALIZE
    catch_hat.initialize();
  #endif
  residuals.allocate(1,nspp,1,nyrs,1,nages,"residuals");
  #ifndef NO_AD_INITIALIZE
    residuals.initialize();
  #endif
  N.allocate(1,nspp,1,nyrs+1,1,nages,"N");
  #ifndef NO_AD_INITIALIZE
    N.initialize();
  #endif
  avail_food.allocate(1,nspp,1,nyrs,1,nages,"avail_food");
  #ifndef NO_AD_INITIALIZE
    avail_food.initialize();
  #endif
  AvgN.allocate(1,nspp,1,nyrs,1,nages,"AvgN");
  #ifndef NO_AD_INITIALIZE
    AvgN.initialize();
  #endif
  F.allocate(1,nspp,1,nyrs,1,nages,"F");
  #ifndef NO_AD_INITIALIZE
    F.initialize();
  #endif
  Z.allocate(1,nspp,1,nyrs,1,nages,"Z");
  #ifndef NO_AD_INITIALIZE
    Z.initialize();
  #endif
  S.allocate(1,nspp,1,nyrs,1,nages,"S");
  #ifndef NO_AD_INITIALIZE
    S.initialize();
  #endif
  M2.allocate(1,nspp,1,nyrs,1,nages,"M2");
  #ifndef NO_AD_INITIALIZE
    M2.initialize();
  #endif
  M2old.allocate(1,nspp,1,nyrs,1,nages,"M2old");
  #ifndef NO_AD_INITIALIZE
    M2old.initialize();
  #endif
  obj_fun.allocate("obj_fun");
}

void model_parameters::userfunction(void)
{
  srv_q = mfexp(log_srv_q);
  M2_pen.initialize();
  M2old.initialize();
  M2.initialize();
  calc_selectivity();
  calc_numbers_at_age(1);
  if (mode)
  {
   for (iter=1;iter<=niter;iter++)
    calc_numbers_at_age(2);
  }
  for (k=1;k<=nspp;k++)
  {
    calc_predicted_values(k);
    calc_objective_function(k);
  }
}

void model_parameters::calc_selectivity(void)
{
  for (k=1;k<=nspp;k++)
  {
    avg_sel(k)               = log(mean(mfexp(sel_coff(k))));
    sel(k)(1,nselages)       = mfexp(sel_coff(k));
    sel(k)(nselages+1,nages) = sel(k,nselages);
    avg_srv_sel(k)            = log(mean(mfexp(srv_sel_coff(k))));
    if (active(srv_sel_slp)) // Asymptotic survey selectivty
    {
      for (i=1;i<=nages;i++)
        srv_sel(k,i) = 1/(1+exp(-srv_sel_slp(k)*(double(i)-srv_sel_inf(k))));  
    }
    else  //survey selectivity uses coefficients
    {
        srv_sel(k)(1,nselages)       = mfexp(srv_sel_coff(k));
        srv_sel(k)(nselages+1,nages) = srv_sel(k,nselages);
        srv_sel(k)                  /= mean(srv_sel(k) );
    }
  }
}

void model_parameters::calc_available_food(void)
{
  for (p=1; p<= nspp; p++)  // predator loop
    {
    for (j=1;j<=nages;j++)  // predator age loop
    {
      for (k=1;k<=nspp;k++)  // prey loop
      {
        for (int prey_age =1;prey_age<=nages;prey_age++)  // prey age loop
        {
          avail_food(p,i,j) += suit_main(p,k,j,prey_age)* AvgN(k,i,prey_age)* wt(p,k,j,prey_age);
       } //end prey age loop
      }  // end prey spp loop
      avail_food(p,i,j) += other_food(p)* suit_other(p,j); 
    } // Pred age (j)
  }  // pred loop
}

void model_parameters::calc_M2(void)
{
  for (int prey_k=1;prey_k<=nspp;prey_k++)  // prey spp loop
  {
    for (int prey_age =1;prey_age<=nages;prey_age++)  // prey age loop
    {
      dvariable Mtmp=0.;
      for (p=1;p<=nspp;p++)   // pred species loop
      {  
        for (j=1;j<=nages;j++)  // Pred age loop
        {
          Mtmp += AvgN(p,i,j) * food_in(p,i,j)*suit_main(p,prey_k,j,prey_age)/avail_food(p,i,j);
        }  // end pred age loop
        M2(prey_k,i,prey_age) = Mtmp;
        //M2(prey_k,i,prey_age)=MIXM2*M2old(prey_k,i,prey_age)+(1.-MIXM2)*M2(prey_k,i,prey_age);
        //M2old(prey_k,i,prey_age)=M2(prey_k,i,prey_age); 
     }  // end pred spp loop
    }   // end prey age loop
  }   // end prey spp loop
 /* for (k=1;k<=nspp;k++)   // pred species loop
  {
    //M2(k,i)    = alpha*elem_prod(AvgN(k,i) , elem_div(pred_food_ration(k,i) , avail_food(k,i)));  //
    M2(k,i)    += elem_div(elem_prod(AvgN(k,i),suit_main( avail_food(k,i)); 
    //M2(K,i) += alpha*AvgN(k,i) * food_in(p,i,j)* suit_transp(p,prey_k,j,prey_age)/ avail_food(p,i,j);
    M2(k,i)    = MIXM2 * M2old(k,i)  + (1.-MIXM2)*M2(k,i);
    M2old(k,i) = M2(k,i);
  }
  */ 
}

void model_parameters::calc_mortality(void)
{
  // if(i!=1)
    F(k,i) =  sel(k) * mean_F(k) * mfexp(F_dev(k,i));
  // else
    // F(k,1) = sel(k)*mean_F(k) * mfexp(F_dev(k,1)) ;
  //cout<<k<<" "<<i<<" "<<M2(k,i)(1,4)<<endl;
  Z(k,i) = F(k,i) + M1(k)(1,nages) + M2(k,i);
  S(k,i) = mfexp(-Z(k,i));
}

void model_parameters::calc_numbers_at_age(int pass_number)
{
  for (k=1;k<=nspp;k++)
     {
     // Sub-vector of rec_dev(k) from 1 to nyrs (rec_dev is longer now...)
     R(k) = mfexp(mean_rec(k) + rec_dev(k)(1,nyrs) );
     // Top row, left column....
     for (i=1;i<=nyrs;i++)
        N(k,i,1) = R(k,i);      // Recruitment  
        // Initial age composition (first year)
        for (j=2;j<=nages;j++)
           N(k,1,j)= mfexp(mean_rec(k) + init_dev(k,j) );
    switch (pass_number)
    {
      case 1:
       AvgN(k,1) = N(k,1);
       break;
      case 2:
      default:
      AvgN(k,1) = elem_div( elem_prod( N(k,1),(1.-S(k,1))) ,Z(k,1));
       break;
     }
  } // End loop over species
  // Main year loop for filling numbers at age
  avail_food.initialize();
  for (i=1;i<nyrs;i++)
     {
      if (mode && current_phase()>2)
        {
         calc_available_food();
         if (pass_number>1)
         calc_M2();
        }
      for (k=1;k<=nspp;k++)
         {
          calc_mortality();
          N(k,i+1)(2,nages)= ++elem_prod(N(k,i)(1,nages-1),S(k,i)(1,nages-1)) ;
          N(k,i+1,nages) += S(k,i,nages)*N(k,i,nages);
          switch (pass_number)
         {
       case 1:
        AvgN(k,i+1) = N(k,i+1);
        break;
       case 2:
       default:
        AvgN(k,i+1) = elem_div( elem_prod( N(k,i+1),(1.-S(k,i+1))) ,Z(k,i+1));
        break;
      }
    }
  }
  i=nyrs;
  if (mode && current_phase()>2 )
  {
    calc_available_food();
    calc_M2();
  }
  for (k=1;k<=nspp;k++)
  {
    calc_mortality();
    N(k,i+1)(2,nages)= ++elem_prod(N(k,i)(1,nages-1),S(k,i)(1,nages-1)) ;
    N(k,i+1,nages) += S(k,i,nages)*N(k,i,nages);
  }
}

void model_parameters::calc_predicted_values(int k)
{
  catch_hat(k) = elem_prod(elem_div(F(k),Z(k)) ,elem_prod(1.-mfexp(-Z(k)) , N(k)));
  // for (i=1;i<=nyrs;i++) for (j=1;j<=nages;j++) catch_hat(k,i,j)=F(k,i,j)/Z(k,i,j) *(1.-mfexp(-Z(k,i,j)))* N(k,i,j);
  for (i=1;i<=nyrs;i++)
  {
    tc_hat(k,i)= sum(catch_hat(k,i));
    p_hat(k,i) = catch_hat(k,i)/tc_hat(k,i);
    biomass(k,i) = N(k,i) * avg_growth(k,i)(1,nages);
  }
  for (i=1;i<=n_srv_biom(k);i++)
  {
    // convert years into indices for 1-=nyrs counting purposes 
    int yr_ind = yrs_srv_biom(k,i) - styr + 1;
    // cout<<k<<" "<<yr_ind<<" "<<nages<<" "<<i<<endl;
    // note: assumes survey occurs in mid-year (pow(,.5))
    srv_biom_hat(k,i) = elem_prod(srv_q(k) * srv_sel(k) , elem_prod(pow(S(k,yr_ind),0.5),N(k,yr_ind)) ) * avg_growth(k,yr_ind)(1,nages);
  }
  for (i=1;i<=n_srv_age(k);i++)
  {
    // convert years into indices for 1-=nyrs counting purposes 
    int yr_ind = yrs_srv_age(k,i) - styr + 1;
    // need to test for type of age data here (if length bins are different
    dvar_vector tmp_age = elem_prod(srv_sel(k) , elem_prod(pow(S(k,yr_ind),0.5),N(k,yr_ind)) );
    if ( srv_age_type(k)==1)
      srv_age_hat(k,i) = tmp_age / sum(tmp_age);
    else
    {
      srv_age_hat(k,i) = tmp_age * srv_age_err(k);
      srv_age_hat(k,i) /= sum(srv_age_hat(k,i));
    }
  }
 // OjO what's up with this?
  if (sd_phase())
  {
    depletion(k) = biomass(k,nyrs)/(1.e-20 + biomass(k,1));
    for (i=1;i<=nyrs;i++)
    {
      // prey_consumed(k,i) = wt(k,i)*elem_prod(M2(k,i),AvgN(k,i)); // Vector * vector = scalar
      //prey_consumed(k,i) = wt(k,i)(1,nages)*elem_prod(M2(k,i),AvgN(k,i)); // Vector * vector = scalar
    }
  }
}

void model_parameters::calc_objective_function(int k)
{
    //residuals(k)= elem_div(obs_catch(k)-catch_hat(k),sqrt(catch_hat(k) + .1));
    // Multinomial
    obj_fun -= tau*sum(elem_prod(p_obs(k),log(p_hat(k)+1.e-4)));
  // Errors in total catch estimation
    obj_fun += 100.*norm2(log(tc_obs(k))-log(tc_hat(k)+1e-4));
    if (!last_phase())
      obj_fun+= 50.*square(log(mean(F(k))/.2));
    else
      obj_fun+= 10.*square(log(mean(F(k))/.2));
  // Invoke a penalty when the partial F's go down with age
    for (j=1;j<=nages-1;j++)
      if (sel(k,j)>sel(k,j+1))
        obj_fun += 20.*square(log(sel(k,j)/sel(k,j+1)));
    obj_fun += 10.*square(avg_srv_sel(k)); // this part will go to zero (it's a condition)
    obj_fun += 10.*square(avg_sel(k)); // this part will go to zero (it's a condition)
    obj_fun += 10. * norm2(first_difference(first_difference(log(srv_sel(k)))));
    obj_fun += 10. * norm2(first_difference(first_difference(log(sel(k)))));
  // cout <<srv_sel<<endl<<srv_q<<endl<<srv_biom_hat<<endl;
  // for (k=1;k<=nspp;k++)
  for (i=1;i<=n_srv_biom(k);i++)
    obj_fun += square(log(srv_biom(k,i)) - log(srv_biom_hat(k,i)) ) / 
                                 (2.*srv_biom_lse(k,i)*srv_biom_lse(k,i));
  for (i=1;i<=n_srv_age(k);i++)
    obj_fun -= srv_age_n(k,i)*sum(elem_prod(srv_age(k,i),log(srv_age_hat(k,i)+1.e-4)));
  obj_fun +=  1.* norm2(rec_dev(k));
  obj_fun +=  1.* norm2(init_dev(k));
  obj_fun +=  1.* norm2(F_dev(k));
  obj_fun += 100.*(square(mean(F_dev(k))) + square(mean(rec_dev(k))) + square(mean(init_dev(k))));
  //if (active(F_dev))
   // obj_fun += 10.* square(trend*F_dev(k))/norm2(F_dev(k)+1e-20);
  //obj_fun += 10.* square(alpha-1);
  obj_fun += 10.* M2_pen;
}

void model_parameters::report()
{
 adstring ad_tmp=initial_params::get_reportfile_name();
  ofstream report((char*)(adprogram_name + ad_tmp));
  if (!report)
  {
    cerr << "error trying to open report file"  << adprogram_name << ".rep";
    return;
  }
  cout <<endl<<"========End of phase: "<<current_phase()<<" ============"<<endl<<endl;
  report << "INPUTS"<<endl;
  report << "cATCH AT AGE"<<endl;
  report << obs_catch <<endl;
  report << "weight at the stomach"<<endl;
  report << wt <<endl;
  report << "food_in"<<endl;
  report << food_in <<endl;
  report << " avg weight"<<endl;
  report << avg_growth <<endl;
  report << "other food"<<endl;
  report << other_food <<endl;
  report << "suit plk as predator, plk as prey, age of prey 1?"<<endl;
  report << suit_main(1) <<endl;
  report << "suit cod as predator, plk as prey, age of prey 1?" << endl;
  report << suit_main(2) <<endl;
  report << "suitabilities other food"<<endl;
  report << suit_other <<endl;
  report << "residual mortality"<<endl;
  report << M1 <<endl;  
  report << "suit transp plk as predator, plk as prey, age of prey 1?"<<endl;
  report << suit_transp(1) <<endl;
  report << "suit transp cod as predator, plk as prey, age of prey 1?" << endl;
  report << suit_transp(2) <<endl;
  //report << "food intake" << endl;
  //report << food_in(1) <<endl;
  //report << "food intake" << endl;
  //report << food_in(2) <<endl;
  report << "OUTPUTS"<<endl;
  report << "Numbers at age" <<endl;
  report << N<<endl<<endl;
  report << "Fishing mortality at age" <<endl;
  report << F<<endl<<endl;
  report << "Survival at age" <<endl;
  report << S<<endl<<endl;
  //report << "Residuals"<<endl;
  //report << residuals  <<endl<<endl;
  report << "Biomass"<<endl<<biomass<<endl;
  report << "Predicted total catch"<<endl;
  report <<  tc_hat <<endl;
  report << "observed total catch"<<endl;
  report <<  tc_obs <<endl;
  report << "Predicted catch at age"<<endl;
  report <<  p_hat <<endl;
  report << "observed catch at age"<<endl;
  report <<  p_obs <<endl;
  report << "Values for M2        "<<endl;
  report << "plk          "<<endl;
  report <<  M2(1) <<endl<<endl;
  report << "pcod          "<<endl;
  report <<  M2(2) <<endl<<endl;
  report << "available food plk"<<endl;
  report << avail_food(1)<< endl;
  report << "available food cod"<<endl;
  report << avail_food(2)<< endl;
  report << "average N plk"<<endl;
  report << AvgN(1)<< endl;
  report << "average N cod"<<endl;
  report << AvgN(2)<< endl;
  report << "======Survey_selectivity-at-age==========="<<endl;
  for (k=1;k<=nspp;k++)
  {
    report << "==Species: "<<k<<"----"<<endl;
    report << srv_sel(k)<<" "<<endl;
  }
  report << "======Survey_Biomass_Fit=================="<<endl;
  for (k=1;k<=nspp;k++)
  {
    report << "==Species: "<<k<<"----"<<endl;
    for (i=1;i<=n_srv_biom(k);i++)
      report << yrs_srv_biom(k,i)<<" "<<srv_biom(k,i)<< " "<<srv_biom_hat(k,i)<<endl;
  }
  report << "======Survey_age_composition_fits=================="<<endl;
  for (k=1;k<=nspp;k++)
  {
    report << "==Species: "<<k<<"----"<<endl;
    for (i=1;i<=n_srv_biom(k);i++)
      report << yrs_srv_age(k,i)<<" "<<srv_age(k,i)<< " "<<srv_age_hat(k,i)<<endl;
  }
  report << "======Fishery_age_composition_fits=================="<<endl;
  for (k=1;k<=nspp;k++)
  {
    report << "==Species: "<<k<<"----"<<endl;
    for (i=1;i<=nyrs;i++)
      report << styr+i-1 <<" "<<p_obs(k,i)<< " "<<p_hat(k,i)<<endl;
  }
  report << "======Total_catch_fits=================="<<endl;
  for (k=1;k<=nspp;k++)
  {
    report << "==Species: "<<k<<"----"<<endl;
    report << tc_obs(k)<< endl<<tc_hat(k)<<endl;
  }
}

void model_parameters::set_runtime(void)
{
  dvector temp("{.1, 0.001,0.00000001}");
  convergence_criteria.allocate(temp.indexmin(),temp.indexmax());
  convergence_criteria=temp;
  dvector temp1("{200, 300, 3000}");
  maximum_function_evaluations.allocate(temp1.indexmin(),temp1.indexmax());
  maximum_function_evaluations=temp1;
}

void model_parameters::preliminary_calculations(void){
  admaster_slave_variable_interface(*this);
}

model_data::~model_data()
{}

model_parameters::~model_parameters()
{}

void model_parameters::final_calcs(void){}

#ifdef _BORLANDC_
  extern unsigned _stklen=10000U;
#endif


#ifdef __ZTC__
  extern unsigned int _stack=10000U;
#endif

  long int arrmblsize=0;

int main(int argc,char * argv[])
{
    ad_set_new_handler();
  ad_exit=&ad_boundf;
  gradient_structure::set_MAX_NVAR_OFFSET(1000);
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(3000000);
  gradient_structure::set_CMPDIF_BUFFER_SIZE(4000000);
    gradient_structure::set_NO_DERIVATIVES();
    gradient_structure::set_YES_SAVE_VARIABLES_VALUES();
  #if defined(__GNUDOS__) || defined(DOS386) || defined(__DPMI32__)  || \
     defined(__MSVC32__)
      if (!arrmblsize) arrmblsize=150000;
  #else
      if (!arrmblsize) arrmblsize=25000;
  #endif
    model_parameters mp(arrmblsize,argc,argv);
    mp.iprint=10;
    mp.preliminary_calculations();
    mp.computations(argc,argv);
    return 0;
}

extern "C"  {
  void ad_boundf(int i)
  {
    /* so we can stop here */
    exit(i);
  }
}
